First-principles molecular-dynamics simulations of a hydrous silica 
melt: Structural properties and hydrogen diffusion mechanism 

Markus Pohlmann 1 ' 2,3 , Magali Benoit 1 and Walter Kob 1 
1 - Laboratoire des Verves, Universite Montpellier II, 34-095 Montpellier, France 
2 - Physik- Department El 3, Technische Universitdt Miinchen, 85747 G 'arching, Germany 
3 - Inistitut Laue-Langevin, 38042 Grenoble, France 

Abstract 

We use ab initio molecular dynamics simulations to study a sample of liquid silica containing 
3.84 wt.% H2O. We find that, for temperatures of 3000 K and 3500 K, water is almost exclusively 
dissolved as hydroxyl groups, the silica network is partially broken and static and dynamical 
properties of the silica network change considerably upon the addition of water. Water molecules 
or free O-H groups occur only at the highest temperature but are not stable and disintegrate rapidly. 
Structural properties of this system are compared to those of pure silica and sodium tetrasilicate 
melts at equivalent temperatures. These comparisons confirm the picture of a partially broken 
tetrahedral network in the hydrous liquid and suggest that the structure of the matrix is as much 
changed by the addition of water than it is by the addition of the same amount (in mole %) of 
sodium oxide. On larger length scales, correlations are qualitatively similar but seem to be more 
pronounced in the hydrous silica liquid. Finally, we study the diffusion mechanisms of the hydrogen 
atoms in the melt. It turns out that HOSi2 triclusters and SiO dangling bonds play a decisive role 
as intermediate states for the hydrogen diffusion. 

PACS numbers: 61.20.-p,61.20.Ja,71.15.Pd,66.10.-x 
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I. INTRODUCTION 



Over the last decades, amorphous silicates containing several wt.% of water have attracted 
the interest of many experimental and theoretical research groups due to its importance 
in science and technology. E.g. these systems are supposed to play an important role in 
magmatic flow in the earth crust and have therefore been investigated extensively by many 
geologists (for reviews see [l|,|2j). In silicate melts (lava), water can be released if temperature 
and pressure conditions are altered and in that case, the released water molecules assemble in 
bubbles. These bubbles can set up pressure for explosive volcanism [0, HI- The microscopic 
origin of such bubbles is still very poorly understood [H and one of the main motivating 
aspects for our study is to find indicators for such bubble formation. The possible use 
of water containing silica in fuel cells [0, |0] and the use of water as defect passivant in 
semiconductor devices (and reladed problems) [7] and its disturbing presence in optic fibers 
[0] underline the technological importance. 

Over the past years pure amorphous silica has been well understood regarding its electronic 
and structural properties and experimental and simulational data seem to agree well [0, HO, 
fill . H2I Fil| . In recent studies these investigations have been extended to the technologically 
highly important class of sodium silicates [3 HO, HO, H3, HO; O New fundamental 
ideas of the structural properties and the diffusion mechanism of sodium atoms in these 
melts have been obtained. 

The impact of water on the structure and dynamics of silica is supposed to be at least as 
important as the addition of sodium [21]. Various methods have been used to investigate 
water speciation and diffusion in a Si02 network. NMR, Infrared and Raman spectroscopy 
as well as neutron scattering experiments suggest a chemical reaction of molecular water to 
SiOH groups in the melt for low initial water concentrations and a saturation of this reaction 
if more water is added to the liquid |l|. m chemical terms it represents a balance of the 
form 

Si - O - Si + H 2 < — ► 2(SiOH) 

that follows the Chatelier principle [22J and which can be shifted to any side by the variation 
of external conditions like temperature or concentration of one species. The energetics is 
supposed to be such that a reaction to the right is endothermal and hence at high tempera- 
tures SiOH is the dominant species p| . On the other hand, the mentioned experiments have 
difficulties to provide reliable quantitative data. The different fabrication processes of the 
samples and the overlap ping of hydroxyl and water signals in the mentioned spectra are the 
main reasons for this |l|, |23[ . In-situ measurements in the geologically relevant liquid state 
are hardly accessible due to elevated temperatures and pressures. 

Computer simulations do usually not encounter these problems since, in principle, they 
can be done at high temperature and pressure. Nevertheless problems for computer sim- 
ulations that use inter atomic potentials arise directly from the above reaction: Potentials 
that describe appropriately all possible types of oxygen (bridging oxygens in the tetrahedral 
network, oxygens in hydroxyl groups, molecular water oxy ge ns ) have been developed with 
some success but do not describe the dynamics very well (2J, |25[ . Here the great advantage 
of ab initio simulations comes into play as the potential and the speciation are found and 
updated from accurate electronic structure calculations at each molecular dynamics time 
step. Especially ab initio Car-Parrinello simulations turned out to be a powerful tool for 
investigations of disordered silica and sodium silicate systems HI HoL fTTl H3 . HO, HO]. For a 
review of previous ab initio simulated amorphous systems see |27j |. 

In this paper we present an equilibrated liquid of water-containing Si02 and its structural 
and dynamical properties. These properties will be compared to that of a silica and a sodium 
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silicate melt when possible. The paper is organized as follows: In section |H] we discuss the 
used simulation techniques and in section ITTT1 we present the structure of the liquid and draw 
some conclusion on diffusion mechanisms for hydrogen. Conclusions are summarized in the 
last section. 



II. SIMULATION DETAILS 



The molecular dynamics simulations were performed with the Car-Parrinello code CPMD 
[iEll^] at 3000 K and 3500 K. The electronic subsystem was treated with a density functional 
(DFT) approach in a generalized gradient approximation using the PBE functional 3(1 IS] • 



The core electrons of silicon and oxygen atoms were treated in a pseudo potential model 
using Troullier-Martins type pseudo potentials [j^. A plane wave T point expansion with 
an energy cutoff of 50 Ry turned out to be sufficient for an appropriate description of the 
inter atomic forces. The 50 Ry cutoff with the PBE functional was first tested on the H2O 
dimer and on a-quartz. These tests showed that e.g. the two SiO distances of a-quartz are 
equal to 1.624 A and to 1.628 A independent of the energy cutoff between 50 Ry and 100 
Ry. On the other hand, we found for the H 2 dimer that the 0-0 distance, the quantity 
which is the most sensitive to a change of the cutoff, shows only a variation from 2.925 A 
to 2.950 A if the cutoff is decreased from 90 Ry to 50 Ry. 

For the Si02-H 2 liquid we used a simulation box containing 102 atoms (30 silicon, 64 
oxygen, 8 hydrogen) corresponding to a concentration of water of 3.84 wt.%. The simulation 
box was prepared by starting from an equilibrated liquid silica configuration (obtained with 
the effective potential proposed by van Beest et al |33[) with 32 Si02 units. Two randomly 
chosen silicon atoms were then substituted by four hydrogen atoms each so that the chemical 
stoichiometry was preserved. These hydrogen atoms were attached to the oxygen atoms 
which were formerly bound to the removed silicon atoms. Using CPMD the liquid was 
subsequently equilibrated at 5000 K for some time and finally cooled down to 3500 K and 
3000 K and annealed for roughly 5 ps before recording the trajectories for the present work. 
Since, to our knowledge, densities of water-containing silica systems (especially in the liquid 
state) have not been measured, the volume of the box had to be chosen such that the 
calculated internal stress of the system had a mean value of zero. This was the case for a 
box size of 11.50 A and hence a density of 2.05 g/cm 3 . 

Comparisons of structural quantities like pair distribution functions, mean square displace- 
ments or angle distributions extracted from molecular dynamics simulations performed at 
3500 K with a cutoff of 50 Ry and of 80 Ry did not show differences within the statistical 
error bars. Thus we concluded that the 50 Ry cutoff was sufficiently high. For equilibration, 
the masses of the ions were all set to 28 a.u. (the mass of a silicon atom). Note that a 
change of the ionic masses does not affect the structure of the liquid since at equilibrium all 
structural quantities are independent of the mass. On the other hand, the increase of the 
ionic masses (from 1 to 28 for hydrogen and from 16 to 28 for oxygen) should increase the 
energy gap in the electronic structure and allow an increase of the Car-Parrinello electronic 
mass and hence the use of a larger time step which thus leads to a faster equilibration. The 
equilibration of the system was performed at the two ionic temperatures of 3000 K and 
3500 K employing Nose-Hoover thermostats and an electronic mass of 600 a.u. (energy x 
time 2 ) at a time step of 4.5 a.u. (0.1088 fs). At high temperature, the electronic gap is too 
small compared to k^T to ensure the decoupling of the ionic and the electronic degrees of 
freedom, which is needed to perform Car-Parrinello dynamics. The use of thermostats is 
therefore compulsory. To speed up the equilibration and to perform an efficient canonical 
sampling, one separate Nose-Hoover thermostat chain for each ionic degree of freedom was 
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used (known as "massive" thermostating |34j). The electrons were controlled with one single 
thermostat chain |35|, |36[. Unfortunately, due to the use of thermostats the direct access to 
dynamical properties is no longer available. 
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FIG. 1: Mean square displacements of the Si, O and H atoms at 3500 K (solid lines) and 3000 K 
(dashed lines). 



The system was equilibrated at two temperatures (3500 K and 3000 K) until in a log-log 
plot the averaged mean square displacements (MSD) of the Si atoms showed at long times 
a slope close to unity. This is illustrated in Fig. ^ where the mean square displacements 
(MSD) of the different atomic species are shown for the 3000 K (dashed lines) and the 3500 
K (bold lines) liquids. Usually MSDs of viscous liquids are composed of three regions: The 
ballistic one in which the atoms move without noticing their neighbors and hence a MSD 
that is proportional to t 2 . This ballistic region is followed by a region where the atoms 
are temporarily confined in a cage made of their nearest neighbors. In this regime, the 
atoms rattle around in the cage without significant displacement, leading to a MSD that 
increases only slowly. Finally the atoms leave this cage and start to show a diffusion motion, 
i.e. a MSD that is proportional to t. The choice of the masses and the thermostats affect 
also the MSD. However, the height of the plateau and the displacement at the onset of the 
diffusional regime should be independent of the thermostat. Hence, we consider the system 
to be equilibrated once the diffusional regime is reached which was the case after 4.4 ps at 
3500 K and 10.9 ps at 3000 K. 

In order to check that the liquids were indeed well equilibrated and that there were no 
aging effects, the trajectories were cut into three equal parts. The averaged mean square 
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displacements were then calculated for each part separately and compared to each other. 
Since the three different averaged MSD did not show any drift, aging effects can be excluded 
and equilibration was indeed obtained after the above mentioned times. All structural 
quantities were extracted from the beginning of the diffusional part up to the total length 
of the trajectory (4.4 to 12.5 ps at 3500 K and 10.9 to 22.5 ps at 3000 K). From Fig. [Owe 
note that at the end of the 3500 K trajectory the MSD is 29.2 A 2 , 8.1 A 2 , and 5.0 A 2 for 
the H, O and Si atoms respectively. At 3000 K the MSDs were 16.4 A 2 , 3.5 A 2 and 1.73 A 2 . 
According to the above mentioned equilibration times and the time step of 4.5 a.u. (0.1088 
fs), the numbers of computed time steps were 114900 at 3500 K and 206800 at 3000 K. 
Using a Hitachi SR8000 computer, one time step takes 56 s on a single processor. 



III. RESULTS AND DISCUSSION 

In this section we present the structural analysis of the liquids at 3000 K and 3500 K. 
Commonly considered quantities like pair distribution functions, bond angle distributions, 
coordination numbers, Q-species distributions and bridging to non-bridging oxygen ratio 
will be discussed and, as far as possible, be compared to the data extracted from ab initio 
simulated silica and sodium silicate (NS4) melts [TBI l37l|. In particular the comparison to 
a recently investigated sodium tetrasilicate melt is highly interesting since the valence 
shell configuration of the sodium and the hydrogen atoms are equivalent. 
In order to analyze the structure with particular attention to the formation and the rupture 
of the silica network, water molecules and hydroxyl groups, it turns out to be useful to 
distinguish several types of oxygen atoms: 

O* hydroxyl group oxygen Si-O-H 

(oxygen with one hydrogen and one silicon nearest neighbor) 
wO water molecule oxygen H-O-H 

(oxygen with at least two hydrogen atoms as nearest neighbors) 
BO bridging oxygen Si-O-Si 

(oxygen with two silicon atoms as nearest neighbors) 
NBO non-bridging oxygen Si-0(-?) 

(oxygen with less than two silicon nearest neighbors) 
03 tricluster oxygen 

(oxygen with three nearest neighbors) 

Si 

03H hydrogen containing tricluster g- >0-H 

(oxygen with two silicon and one hydrogen neighbor) 

Note that, within these definitions, all the O* are also counted as NBO and that all the 03H 
are also counted as 03 as well as BO. The wO oxygen atoms can also have silicon neighbors, 
e.g. the oxygen atom in H 2 OSi is considered as a wO. 

The Si and H nearest neighbors of the oxygen atoms were determined as the Si and H atoms 
being located within a sphere the radius of which is given by the positions of the first minima 
in the Si-0 and O-H pair distribution functions, respectively (see section IHI A|) . 



5 



A. Radial Distribution Functions 



In this subsection we analyze the short range correlations of the liquid in terms of the radial 
distribution functions (RDF) 



D aP {r) 
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and A/3 = N/3 — 5 a 



(1) 



where D a p(r) is the number of inter atomic distances between a and (3 atoms found between 
r and r + dr. N denotes the number of atoms of each kind i and V is the volume of the 
simulation box. The corresponding integrated coordination numbers (ICN) are 



(2) 



ICN Q/3 (r) = — / D aP (r') dr' with A a = N a - 6 a/3 
A n Jo 



The RDF and ICN are presented in the left panel of Fig. |21 for the network forming atoms 
silicon and oxygen (a, (3 =Si,0) and in the right panel for the pairs involving H. 
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FIG. 2: Radial distribution function (solid lines) and integrated coordination numbers (dashed 
lines) of Si-Si (a), Si-0 (b), 0-0 (c), Si-H (d), H-0 (e), and H-H (f) for the Si0 2 -H 2 liquids at 
3500 K (bold lines) and 3000 K (thin lines). 

From this figure, one recognizes that the distributions become broader as the temperature 
is increased. The RDFs involving the matrix atoms show a first peak followed by a well 
defined first minimum which becomes more pronounced if the temperature is decreased 
from 3500 K to 3000 K. In particular, the Si-0 RDFs present, after the first peak, a very 
well defined minimum at 2.37 ± 0.05 A for 3500 K and at 2.35 ± 0.05 A for 3000 K. From 
the position of the Si-0 first peak, we can deduce that the most probable Si-0 distance is 
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around 1.65 ± 0.02 A for 3500 K and around 1.63 ± 0.02 A for 3000 K. 
The ICN for Si-0 exhibits a plateau at a value of 4 which indicates that every silicon 
atom has on average four oxygen neighbors and hence that the principal units - the SiC>4 
tetrahedra - are preserved also in the presence of water. Indeed only a small percentage of 
threefold and fivefold coordinated Si atoms are found in the liquids (7 % fivefold corrdinated 
and 4 % threefold coordinated at 3500 K and 2 % fivefold coordinated and 1 % threefold 
coordinated at 3000 K). On the other hand, the ICN for Si-Si shows an inflection point at 
around 3.6, a value which is smaller than the one for a perfect tetrahedral network, 4.0, 
indicating that the tetrahedral network is partially broken. 

Comparing the matrix distributions to those of the pure silica melt, we note that the addition 
of water does not alter signicantly the shapes of the Si-O, 0-0 and Si-Si RDFs presented 
in (3?J , as it was already the case in the sodium silicate melt [l5[ upon the addition of sodium. 

The RDFs for H-Si and H-O, right panel of Fig. |2l show somewhat better defined inter 
atomic distances as the temperature is decreased, as it was the case for the Si-O, Si-Si and 
0-0 RDFs. On the other hand the H-H distribution seems to deviate from this behavior. 
Concerning the H-Si distribution, the first maximum is found around 2.3 A. Since this value 
is much larger than the most probable Si-O distances (1.65 ± 0.02 A and 1.63 ± 0.02 A ), 
the presence of stable molecular Si-H units is excluded. The height of the first peak in Si-H 
and the absence of a well-defined minimum reflect the presence of a broad distribution of 
these distances in the liquid. In contrast, the well-defined first peak in the H-0 distribution 
functions, located at 0.99 ± 0.01 A, followed by a well defined minimum at 1.48 ± 0.03 A 
at 3500 K and at 1.43 ± 0.03 A at 3000 K, reveals the O-H bond as the dominant stable 
configuration for the hydrogen atoms. 

Because of a lack of a well-defined first peak in the RDF for H-H, the existence of stable H2 
molecules can be excluded too. The striking point in the H-H distribution is the difference 
between the RDF at 3000 K and 3500 K between 1.5 A and 3.0 A. Whereas at 3000 K 
there are almost no H-H pairs with a distance around 2.0 A we find a weak peak at this 
distance at 3500 K. This distance is very close to the H-H distance in a free water molecule 
and hence we have evidence that at this temperature water molecules do exist. Indeed, in 
a more detailed analysis of the coordination numbers as a function of time, we found that 
oxygen atoms with two nearest neighbor hydrogens were stable over times of the order of 
500 fs. The absence of this phenomenon at 3000 K is probably due to the considerably 
lower displacement of the hydrogen atoms (see Fig. ^) at this temperature which, because 
of a too small trajectory length, prevents to find two hydrogen atoms sufficiently close to 
each other to form a water molecule. 

We recall at this point that the configuration of the valence shell is the same for the 
hydrogen and the sodium atoms (although it is well known that the physical and chemical 
properties of compounds of equivalent stoichiometry involving these two atom types (such 
as H 2 and Na20) are rather different). This behavior motivates a comparison of the 
results of this study to those of the liquid sodium silicate [HI]. In Fig. El we present 
the comparison of the above discussed hydrogen-containing RDFs with the corresponding 
Na-containing RDFs in the sodium tetrasilicate melt at 3500 K from Ref. [15j]. Two main 
properties seem to govern the character of the distributions: The size of the atoms and the 
ability of the two atom types to form either covalent or ionic bonds. In particular, the RDF 
for Si-X and X-X (X = Na, H) shown in Fig. |3] look quantitatively quite similar except 
that the peaks are shifted to larger distances in the case of Na. This behavior can be easily 
related to the atom size. In contrast to this, the O-X distribution function for X=H has a 
very different shape from the one for X=Na. As described above, for the O-H correlation, 
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FIG. 3: Si-X, X-0, X-X (X=H,Na) radial distribution functions (solid lines) and integrated 
coordination numbers (dashed lines) for the hydrous silica (bold lines) and the NS4 melt from Ref. 
Q (thin lines) at 3500 K. 



a sharp peak at an O-H distance of 0.99 A ± 0.01 is followed by a well-defined minimum. 
In contrast to this, the first peak of the O-Na RDFs is much broader with a maximum at 
2.32 ± 0.05 A followed by a shallow minimum around 3.6 A. These latter differences are 
the signature of the strong covalency of the OH bond and the ionicity of the Na-0 bond. 

To conclude the section on the radial distribution functions, we want to investigate how the 
different oxygen types defined above contribute to the RDF. In order to count the number 
of Si and H neighbors of a given oxygen atom, we used cutoff distances extracted from the 
first minima of the Si-0 and H-0 radial distribution functions, respectively (cutoff for SiO 
= 2.37 A and 2.35 A and cutoff for HO=1.48 A and 1.42 A for the higher and the lower 
temperature, respectively). Figure 0] shows the H-0 radial distribution functions for the 
different oxygen types (O*, NBO, BO) compared to the total H-0 RDF at 3500 K. Results 
obtained for 3000 K are very similar and are therefore not shown here. 

We note that the RDFs for all the different oxygen types show a pronounced peak around 
1.0 A but that the height of this first peak depends strongly on the considered oxygen type. 
The H-O* RDF exhibits a first peak of height 65 which is much larger than the one for the 
total H-0 which is around 8. In contrast to this, the RDF for H-BO has a first peak of height 
1. The most striking point in Fig. 0]is that the peak for H-O* is significantly higher than the 
one for H-NBO. Since all O* atoms are also NBO atoms (remember that by definition the 
O* atoms have a hydrogen atom as second neighbor whereas for the NBO atoms this is not 
necessarily the case), the presence of a large number of Si-0 dangling bonds becomes now 
evident. We will discuss the existence and the temperature dependence of these dangling 
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bonds below. From the inset of Fig. it can be concluded that the H-NBO contribution 
is dominant for distances between 1.5 A - 2.5 A. This means that the NBO atoms are also 
connected to hydrogen atoms via the so-called hydrogen bonds the length of which is equal 
to 2.0 A in the water dimer. Though also the RDF for the other oxygen types show a 
contribution in this range, it is less important. 

In Fig. we present the contribution of the different oxygen types to the Si-0 radial 
distribution function. The height of the first peak of the Si-BO RDF is around 9.3 whereas 
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FIG. 5: Pari distribution functions at 3500 K for Si-0 (bold solid line), Si-BO (dotted line), 
Si-NBO (dashed line), and Si-O* (thin solid line), (b) and (c) are magnifications of the first peaks 
and first minima, respectively. 



the peak height of the total Si-0 RDF is close to 8.7 which indicates that the correlation 
of the bridging oxygen with silicon atoms is stronger than the total Si-0 correlation (see 
Fig. EJd). The heights of the first peak in the Si-O* and Si-NBO RDFs are very close to 
each other at a value of ~ 5.1. As shown in Fig. EJo, a slight shift of the peak positions 
can be observed: The Si-NBO and Si-O* peak positions are around 1.60 A and 1.62 A, 
respectively, whereas the Si-BO peak position is located at 1.63 A. This result implies 
that the tetrahedra having one or more NBO atoms are distorted, as it has already been 
observed in other silicate systems (see Ref. 3 and references therein). 
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Figure |U shows the presence of a considerable jump in the RDFs for Si-O* and Si-NBO 
at the Si-0 cutoff. Between 2.37 A (the Si-0 cutoff) and 4.0 A, the Si-O* and Si-NBO 
RDFs are larger than the total one (Fig. [HJj). This jump and the dominance of the Si-O* 
and Si-NBO can be associated to oxygens of the O* type or to NBOs that are connected 
to a second silicon atom by a weak Si-0 bond. As soon as the distance to the second 
silicon atom becomes smaller than the cutoff (i.e. an 03H or BO is formed) the oxygen 
is considered as BO and therefore the Si-O* or Si-NBO distribution drops abruptly to 
zero at the cutoff distance. Thus this behavior gives insight into the formation and decay 
processes of the 03H clusters and the Si-0 dangling bonds, that will turn out to be one of 
the essential transition states for hydrogen diffusion. 



B. Angular Distributions 



Figure |H1 presents the different angular distributions in the Si02-H20 liquids at 3500 K and 
3000 K and the comparison with NS4 at 3500 K. We remark that the angular distributions 
for the network formers (Si-O-Si) and (O-Si-O) are close to the ones of pure silica at 3500 
K (3?]]. In agreement with the behavior of the RDFs (Sec. 1111 A|) . the angular distributions 
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FIG. 6: (a)-(c) Angular distributions of the hydrous silica sample at 3500 K (bold lines) and 3000 
K (thin lines); (d) comparison to NS4 at 3500K (dashed line). 



at the lower temperature are sharper and more peaked. The Si-O-Si distribution shows 
an additional hump at 90° that is also present in pure silica at 3500 K but much less 
pronounced. This hump is due to the two membered rings which are formed by two 
Si04 units connected by an edge and having two common oxygen atoms. As will be 
seen in Sec. IIII Dl these units are related to intermediate states that are relevant for the 
diffusion process. From the Si-O-Si distribution in Fig. |Hl we note that the number of such 
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rings increases with increasing temperature. The O-Si-0 angular distributions present a 
well-defined maximum at « 110° which corresponds to a typical intra-tetrahedra angular 
distribution (ideal tetrahedra angle ~ 109°). Finally, we can recognize that the temperature 
dependence of the Si-O-H distribution seems to be significantly weaker than the one of the 
Si-O-Si and O-Si-0 distributions. 

The Si-O-X (X= H, Na) angular distributions at 3500 K (Fig. 01) present very similar 
shapes with a slight shift to larger angles for the sodium silicate liquid. In both cases we see 
a broad maximum at angles between 108° and 115°. A maximum in this range shows that, 
in both cases, the chemistry of a twofold coordinated oxygen with a lone pair is realized. 
The shift of the mean angle to a higher value for the sodium silicate is certainly an effect of 
the larger size of the sodium atom since for larger atoms the repulsion of the electron shells 
of the two oxygen-ligand atoms becomes important. 



C. Structure Factors 

We calculated the neutron scattering structure factor according to 
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where b^ (fc=Si,0,H) are the coherent neutron scatte ring lengths and (•) is the thermal 
average. The scattering lengths were taken from Ref. (38( where b^ = 0.41491 x 10~ 14 m, 
bo = 0.5803 x 10~ 14 m and &h — —0.374 x 10 -14 m are reported. A very important feature for 
the experimental verification of the simulation is the replacement of hydrogen by its isotope 
deuterium. Since the coherent scattering length of deuterium (do — 0.671 x 10~ 14 m) is 
quite different from that of hydrogen, the contribution of the hydrogen atoms to the total 
structure factor can be revealed. Figure [?K shows the neutron scattering structure factor 
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FIG. 7: (a) Total neutron scattering structure factor of the hydrous silica liquids at 3000 K (thin 
line) and at 3500 K (bold line) and of a pure silica liquid at 3500 K (dashed line); (b) Total neutron 
scattering structure factor for the hydrous silica liquid at 3500 K (bold solid line), the deuterated 
liquid (bold dashed line) and the difference (thin solid line). 



S n (q) for both simulated temperatures in comparison to the simulated dry silica liquid at 
3500 K [37]. For q > 2.0 A" 1 , the three curves can be considered to be identical which is 
in agreement with the fact that we did not see any differences in the Si-Si, Si-0 and 0-0 
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radial distribution functions for silica, the sodium silicate, and the hydrous silica liquids. 
In contrast to this, the q < 2.0 A _1 region of S n (q) of the hydrous silica liquid exhibits a 
prepeak at 0.95 A" 1 at 3500 K and at 1.3 A -1 at 3000 K which is not present in the pure silica 
melt. (We mention that the box size of the hydrous silica liquid is 11.5 A, corresponding 
to a minimum value of q of 0.55 A -1 and thus the prepeak found at 0.95 A -1 can not be 
attributed to a size effect. Note also that the number of silicon and oxygen atoms was 
almost equivalent in the simulation of pure silica and hence the statistical accuracy is as 
well comparable.) 

For the sodium silicate it turned out that the long range correlations (q < 2.0 A -1 ) contain 
hig hly important information on the diffusion mechanism as it was shown recently 0, 
llTL l39j |. For the sodium silicate it turned out that the long range correlations (q < 2.0 
A -1 ) co ntain important information on the diffusion mechanism as it was shown recently 
[iH Fnl l39l | . It was found that in the liquid a network of channels is formed which enables 
the sodium atoms to move rapidly through the SiO-matrix. This channel structure is thus 
able to explain the relatively high diffusion constant of sodium atoms in sodium silicate 
melts. The characteristic distance between these channels is around 6 A, i.e. two tetrahedra 
diameters, and gives rise to a prepeak in the structure factor at about 1 A" 1 , a structural 
feature which has indeed be found in recent neutron scattering experiments [39]. 
Fig. [7)d shows the structure factor at 3500 K for the hydrated silica and deuterated silica 
liquids as well as their difference. We note remarkable differences of the two signals at 0.95 
A" 1 , around 2.0 A -1 and around 5.0 A -1 . Since the most important contributions to the 
difference are located in the g-vector region of the prepeak (at 0.95 A -1 ), the prepeak can 
be directly attributed to the presence of the hydrogen atoms. For a deeper understanding of 
the prepeak origin in the hydrous melts, we analyzed the partial structure factors presented 
in Fig. El 

The partial structure factors S a p(q) are related to the total S n (q) presented in Fig. [7| as 
follows: 

1 N 

S n (q) = h2 Y, b ^h s a /3(q) (4) 

where S a p(q) are given by 

^) = ^££<exp[iq-(r fc -r,)]} (5) 
iV k 1 

and f al 3 is equal to 0.5 for a 7^ (3 and equal to 1.0 for a = (3. The partial structure 
factors in Fig. |H] can be separated in a group describing the silica matrix distributions 
(Si-Si, Si-O, O-O) and in a group involving distributions with H (Si-H, O-H, H-H) as 
we did for the RDFs (see Sec. IIII A|) . For the matrix part, we find for q > 2.0 A _1 
a nearly perfect agreement between the structure factors presented in Fig. |S] and the 
ones extracted from the NS4 liquid simulation from Ref. ^5[. The only difference is 
the prepeak at 0.95 A -1 for 3500 K and at 1.3 A -1 for 3000 K. We recognize therefore 
a modification of the matrix at a length scale between 4.8 A and 6.6 A which was not 
visible in the RDF representation. Since the prepeak is present not only in the partial 
structure factors involving hydrogen, but also in the 0-0 and Si-Si distributions, we 
note that the network modification does concern not only the hydrogen atoms. This 
behavior can be understood by taking into account the strong and well defined bonding of 
the hydrogen atoms to the silica network (see the O-H radial distribution function in Fig. |2J). 
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D. Distribution of Nearest Neighbors 

In this section we discuss the nearest neighbor coordination numbers. The nearest neighbors 
of an atom were again determined as the atoms that are located within a sphere with a radius 
that is given by the positions of the first minima in the corresponding radial distribution 
function. Whereas the average coordination number has already been extracted from the 
RDFs themselves (Fig. |2J), we present now the distributions. The discussion is limited to 
the coordinations of the network forming atoms. The errors for the results presented in this 
section are related to the number of independent configurations in the liquid. We assume a 
configuration to be independent of a previous one if they are separated by the time it takes 
for the atoms to show a diffusive motion. From Fig. ^we thus recognize that for T=3000 
K and T=3500 K we have two and three independent configurations, repectively. Thus the 
relative errors are l/y/2 = 0.71 at 3000 K and l/y/3 = 0.58 at 3500 K and are then divided 
by the square root to the considered average quantities. The resulting absolute errors will 
be given in the text. 

As for the other quantities described above, we begin with the temperature dependence of 
the distributions in the hydrous liquid and extend the study of the comparison of these 
distributions with the ones found in the sodium silicate liquid. 

Figures Ok and b present the Si-0 and O-Si distributions for the two temperatures of the 
hydrous silica liquid. Both coordination probabilities confirm the picture drawn from the 
analysis of the previous quantities: The silica network is still present in the hydrous silica 
liquid since we find a maximum of the Si-0 distribution at a coordination of four and a 
maximum of the O-Si distribution at a coordination of two. 
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FIG. 9: Si-0 (a) and O-Si (b) coordination distributions for the hydrous silica liquids at 3500 K 
(black bar) and 3000 K (white bar) and Si-0 (c) and O-Si (d) for the sodium tetra-silicate at 3500 
K (grey bar). The black bars in (c) and (d) are identical to the ones in (a) and (b) and represent 
the hydrous liquid at 3500 K. 

In agreement with the temperature dependence of the RDFs and of the angle distributions, 
the coordination distributions become broader as the temperature is increased. For the 
higher temperature, the Si-0 distribution shows significant contributions at coordination 
numbers of three and five and the O-Si distribution shows contributions at a coordination 
of one. 

In addition, the O-Si coordination reveals the speciation of the water molecules, since we do 
not find a significant contribution at a coordination number of zero. Therefore the existence 
of free water molecules can be excluded (a free water molecule has zero silicon neighbor). 
Note that in reality the numerical value is non-zero (but too small to be visible in the figure) 
because of the presence of transient water molecules in the liquid at 3500 K. In contrast, 
the contribution at a coordination of one is roughly | at both temperatures, corresponding 
to the fraction of the non-bridging oxygen in the system that compensate the charges of the 
eight hydrogen atoms in the form of O-H groups. 

Since the dissolution product of water in silica is almost exclusively made of Si-OH units, the 
dissolution mechanism is similar to that of disodium oxide in the NS4 liquid. To compare 
these two mechanisms, we present the Si-0 and O-Si coordination numbers of the hydrous 
silica liquid and of the sodium tetrasilicate liquid at 3500 K in Figs. 0t and d. When 
comparing these two systems, the concentration of sodium atoms of 13.3 mol % in the 
NS4 melt and the concentration of hydrogen atoms of 7.8 mol % in the hydrous silica melt 
should be taken into account. The Si-0 coordination distribution has basically the same 
shape for both systems with a clear maximum at a coordination number of 4 even though 
the distribution for NS4 is, however, somewhat broader. The O-Si coordination probability 
presents stronger differences. Both distributions have their maximum at 2 (indicating that 
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in both cases the silica network is still present) but the absolute values are quite different. 
In both cases the O-Si probability for having one silicon neighbor deviates significantly from 
zero. This indicates the formation of O* in the hydrous melt and equivalent oxygen types 
in the NS4 melt (hereafter denoted by NaO*). If the dissolution product was exclusively 
given by O* and NaO*, the coordination probability would correspond to the Nh/No ratio 
or to the NN a /No ratio, respectively, where Nh, No and Nn & are the number of hydrogen, 
oxygen or sodium atoms in the different systems. These ratio are equal to || = 0.222 for 
the sodium silicate liquid and to ^ = 0.125 in the hydrous silica one. Figure Eli shows a 
contribution of 0.139 ± 0.027 for the hydrous silica sample and of 0.218 ± 0.037 for the NS4 
sample which is in agreement with the simple theoretical prediction. In the hydrous case, 
we note again the important presence of Si-0 dangling bonds at 3500 K compared to the 
NS4 liquid at the same temperature. 

A quantity closely related to the O-Si and Si-0 coordinations is the distribution of the Q n 
species, where n denotes the number of bridging oxygens attached to a silicon atom. For the 
simulated systems, we note that for a perfect dissolution of the water into OH groups, one 
expects a probability for Q 3 equal to = 0.267, the ratio between the number of hydrogen 
and silicon atoms, if all OH groups are attached to different silicon atoms. Indeed Fig. El 
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FIG. 10: Probability to have a Q n speciation for the hydrous sample at 3500 K (black bars) and 
3000 K (white bars). 

exhibits a probability of Q 3 sites of 0.229 ± 0.062 at 3000 K and of 0.269 ± 0.055 at 3500 K, 
respectively. The contributions at n = 2 of 0.029 ± 0.022 at 3000 K and of 0.023 ± 0.016 at 
3500 K are relatively small. In particular, these contributions are smaller than Jj^j = 0.062, 
the probability that two H among the eight ones of the system are found on the same Si 
tetrahedra, which indicates that the Si-O-H groups tend to avoid each other. Due to the 
absence of contributions at n — 1 and n = 0, we exclude the possibility for a clustering of 
O-H groups on specific silicon atoms. Consequently a relatively homogeneous distribution 
of the O-H groups over the silicon atoms is observed. Besides, note that one Q 2 site in the 
system would give a probability of = 0.033 which is larger than the probabilities for Q 2 
presented in Fig. El This indicates that the Q 2 sites do not exist all along the trajectories. 
The Q n species distribution is experimentally accessible with NMR spectroscopy. Farnan 
et al. have determined these quantities for hydrous silica glasses from 29 Si NMR spec- 
troscopy. Their samples contained 2.5 and 8.7 wt.% H2O and the measured Q 3 probabilities 
were 0.086 ± 0.008 and 0.235 ± 0.025, respectively, e.g. values that are significantly lower 
than the Q 3 probabilities found in the present study. This might be due to the fact that 
in the samples of Ref. [40] only half or less of the water molecules were dissolved into OH 
groups. For the 2.5 wt.% H2O sample, the authors found a ratio of OH over H2O of 1/1 
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and an even lower value in the 8.7 wt.% H 2 sample. Since the samples were prepared at 
only 1550°C, these difference in the concentration of the Q 3 species might be due to the 
difference between the experimental temperature and the present simulation ones. 
We also evaluated the probabilities to find oxygen atoms of different types during the tra- 
jectories at 3000 K and 3500 K. In particular the appearance of intermediate states like Si-0 
dangling bonds and 03H triclusters as well as water oxygens (wO) can have consequences 
for the hydrogen diffusion in the melt. Figure fTTI shows the probability to find O* and 03H 
cluster units in the trajectories at 3000 K and 3500 K. 
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FIG. 11: Probability to find O* (a) and 03H cluster (b) in the hydrous silica liquid at 3500 K 
(black bar) and 3000 K (white bar). 



We recall that a dissolution of all hydrogen atoms into O-H groups would give an exact 
number of eight O* (the number of H atoms). As it can be seen in Fig. ITTk . the probability 
distribution is shifted towards values smaller than eight and that at 3500 K the maximum is 
not even at eight. In contrast, the number of 03H cluster increases with increasing temper- 
ature (Fig. ITTb). Since we find that the number of free O-H groups and free water molecules 
is negligible, we conclude that the decrease in concentration of O* at high temperature is 
due to a direct conversion into 03H clusters. The important re-decay of 03H clusters into 
O*, which will be presented in Sec. IIIIE1 seems to confirm this hypothesis. From Fig. Eland 
the simulation of pure silica at 3500 K |37|. we know that three silicon coordinated oxygen 
atoms are hardly present. The tricluster site formation is hence facilitated by the presence 
of hydrogen atoms. Obviously the O-H groups constitute "dead-end-pieces" and a higher 
angular mobility for an O* oxygen than for a BO one can be expected. This higher angular 
mobility enables the approach of the O* to another silicon atom and hence the formation of 
a tricluster. The lower mobility of BO atoms seems to suppress this process in pure silica. 
We have also studied the relation of the 03H cluster and the two membered rings and found 
a significant correlation since more than half of the 03H clusters is part of a two membered 
ring. 

The existence of Si-0 dangling bonds has been demonstrated in Fig. 0] In Fig. ^]we present 
the temperature dependence of the probability to find a given number of these oxygen types 
in the hydrous silica liquid. 

Whereas at 3000 K the maximum of the distribution is at zero, it shifts to a value of two 
at 3500 K and the probability to find zero dangling bonds at this temperature is small. 
Hence, the formation of these species is an effect of the elevated temperature. Since this 
structural feature is not found in a pure silica liquid, the question of its origin emerges. We 
relate the existence of these dangling bonds directly to the formation of the 03H cluster. 
Since these tricluster units increase the coordination of the oxygen atom from two to three, 
the system attempts to compensate of this stoichiometry violation. Since the coordinations 
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FIG. 12: Probability distribution of the Si-0 dangling bonds for the hydrous silica liquid at 3500 
K (black bars) and 3000 K (white bars). 



of the hydrogen and silicon atoms are rarely violated, the appearance of the Si-0 dangling 
bonds (where O has only one neighbor) is certainly associated to the appearance of the 03H 
clusters. 

Now that the existence of the Si-0 dangling bonds is evident, the question of their contri- 
bution to an eventual fast hydrogen diffusion emerges. Do these dangling bonds serve as 
acceptor and donator states for hydrogen ? From Figs. H]and|Hl we know that Si-0 dangling 
bonds exist both with a weak bond to a silicon atom and with a weak bond to a hydrogen 
atom. But the analysis of the recombination of these sites shows that only about 17 % of 
the dangling bonds recombine to a O* site and 83 % recombine to a BO site (see Sec. 1111 Ej) . 
The last oxygen species we want to present here is the water oxygen (wO). As already 
discussed in the introduction, these units are supposed to play a decisive role for the proton 
transport. We also know from the H-H RDF first peak at 3500 K (Fig. EJ) that such units 
exist at this temperature. Figure [TBI shows the probability to find such water units. Indeed 




number of H 2 (wO) units 

FIG. 13: Probability for H 2 units at 3500K (black bars) and 3000K (white bars). 

at the higher temperature, one water unit is present with a probability of roughly 10 % 
along the trajectory. At 3000 K this probability is only 2 %. The corresponding water 
concentrations are 0.30 mol % and 0.06 mol %, respectively. Note that the concentration of 
0.30 mol % at 3500 K should give rise to a probability of 0.003 for the zero O-Si coordination 
in Fig. El if the water molecules were free. Since the contribution at zero is two orders of 
magnitude lower, we underline the important point that those water units are not free. A 
detailed analysis reveals an 03 coordination with two hydrogen neighbors and one silicon 
neighbor. As already mentioned above, the found water unit concentrations are far remote 
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from experimentally measured water concentrations. Considering again the data of Farnan 
et al. |4£| where a ratio of OH over H 2 of 1/1 (50 mol % H 2 0) at a total water content of 
2.5 wt.% is found, one would expect important contributions around two in Fig. EH Again 
these differences in the concentration of molecular water probably arise from the difference 
between the experimental temperature and the present simulation ones. 

We conclude the section on the structural properties by a brief summary of the fig- 
ures that showed evidence of the existence of the different intermediate states and that 
quantify their probability of occurence in the liquid: 



structural unit 


relevant figures 


SiOH groups 
Si-0 dangling bonds 
03H triclusters 
water molecules 


2e (evidence), 11a (quantification) 

4 (evidence), 5 (evidence), 12 (quantification) 

5c (evidence), lib (quantification) 

2f (evidence), 13 (quantification) 



TABLE I: Figures that show evidence for the existence of the transition states and that quantify 
their probability of occurence. 



E. Hydrogen Diffusion Process 

In this section we discuss the possible mechanisms for hydrogen diffusion as it has been 
mentioned in the introduction. We emphasize again that, due to the presence of thermostats, 
dynamical quantities like the diffusion constants cannot be extracted reliably. Nevertheless, 
the structural characteristics of the melt should allow to obtain at least some insight into 
its dynamical properties. Our aim is therefore to determine whether the structural units we 
discussed in Sec. IIII D[ such as the 03H clusters or the Si-0 dangling bonds, can serve as 
intermediate states for the hydrogen diffusion processes in Si0 2 -H 2 liquids. 
We start this discussion by making a list of the possible hydrogen diffusion mechanisms 
in liquid silica, eliminating the free water molecules or stable O-H groups as possible free 
hydrogen carriers since, as discussed above, they are absent in our simulation data. As the 
hydrogen atoms are attached to the silica matrix in the form of Si-O-H groups, three possible 
mechanisms come into play: 

1. motion of the hydrogen in the form O-H O — > O H-0 

2. motion of the oxygen in the form H-0 H — ► H O-H 

3. motion of the O-H group in the form Si-(OH) Si — > Si (OH)-Si. 

The first two mechanisms require the rupture of an O-H bond whereas the third one requires 
the rupture of a Si-0 bond. A typical reaction involving the two first processes are shown 
in Figs. HU and [T3] . 

In the simulated trajectories, we can easily count the number of O-H ruptures in order to find 
out whether the first two mechanisms exist in the liquid. However, it is possible to distingish 
the two processes only by looking at the reactants and the decay products associated to each 
encountered O-H rupture. The possible reactants (resp. decay products) for the formation 
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FIG. 14: Typical hydrogen diffusion reaction of process 1.. A hydrogen is released from an Si-O-H 
(O*) to a bridging oxygen (BO) forming subsequently another O*. A Si-0 dangling bond and an 
unsaturated silicon atom are the resulting products. 




FIG. 15: Typical hydrogen diffusion reaction of process 2.. A hydrogen is released from an Si-O-H 
(O*) to another Si-O-H group (O*) forming a water like structure. An Si-0 dangling bond and an 
relaesed hydrogen atom are the resulting products. 



(resp. rupture) of an O-H group are given in Tab. ITT1 where the O* refers to an oxygen in 
a Si-O-H group, wO to an oxygen atom in a water molecule, BO to an oxygen atom in a 
Si-O-Si group and 03H to an oxgen atom in a Si-(OH)-Si group (see the definitions in the 
introduction of Sec. IHIjl . Note that the formation or rupture of an O-H bond associated 
with a wO — > O* requires the water molecule to be close to a silicon atom which is indeed 
the case since we have deduced from Fig. El that no free water molecule exist in the liquid. 



O-H formation 


O-H rupture 


reactant final state 


initial state decay product 


Si-0 dangling 0* 


0* Si-0 dangling 


0* wO 


wO 0* 


BO 03H 


03H BO 



TABLE II: List of possible reactants for the O-H formation and of possible decay products for 
the O-H rupture. The notations refer to the definition of the different oxygen types given in the 
introduction of Sec. II I II The final (resp. initial) state gives the type of the oxygen atom in the 
formed (resp. destroyed) O-H unit. 

The average number of O-H bonds in the Si02-H20 melts is found to be equal to 8.4 at 3500 
K and to 8.1 at 3000 K and is therefore larger than the total number of hydrogen atoms in 
the system. Hence we conclude that intermediate states with two oxygen atoms close to a 
single hydrogen exist and that they must serve as intermediate states for hydrogen exchange 
reactions (process number 1.). 

By counting the number of O-H ruptures along the trajectories and the type of reac- 
tants/decay products associated to these ruptures, we are able to show the existence of 
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hydrogen diffusion processes of types 1. and 2.. For the following we assumed that an 
O-H bond was formed if the O-H interatomic distance was larger than the O-H cutoff at 
step t and smaller than the O-H cutoff at step t + At (At being the length of a time step 
of 0.1088 fs). The nature of the reactant is found by counting, at time t, all silicon and 
hydrogen neighbors (again in terms of the nearest neighbor cutoff) of the later O-H oxygen. 
A rupture of an O-H bond was assumed to have happened if one H and one O atom had an 
interatomic distance smaller than the O-H cutoff at time t and if their interatomic distance 
became larger than the O-H cutoff at step t + At, irrespective of whether subsequently the 
same bond was formed again or not. The decay product of an O-H rupture is the H donator 
unit (without the H atom) and is found by counting, at step t + At, all silicon and hydrogen 
neighbors (again in terms of the nearest neighbor cutoff) of the former O-H oxygen. 
At 3500 K the total number of O-H formations and ruptures was 97 and 94 respectively, 
during the equilibrated part of the trajectory (8.1 ps). At 3000 K we found 68 O-H for- 
mations and 69 ruptures during the equilibrated part of the trajectory (11.6 ps). However, 
in order to take into account the O-H formations and ruptures that serve for the hydrogen 
diffusion, we also have to distinguish between the formations and ruptures that give rise to 
hydrogen transfers or recombinations. In order to perform this separation between transfers 
and recombinations it is necessary to relate each decay of an OH bond to a formation of an 
OH bond. A hydrogen transfer requires the rupture of one O-H bond and the formation of 
a different one, whereas a recombination implies only one bond. Note that the rupture of 
the first bond will occcur before the formation of the new bond if a transient free hydrogen 
is formed, or, as in almost every case, after the formation of a new bond which implies the 
formation of the intermediate state. Making this distinction between transfers and recombi- 
nations, we find 28 transfers at 3500 K (corresponding to a ratio of || = 28.8 % transfers and 
71.2 % recombinations) and 8 transfers at 3000 K (corresponding to a ratio of = 11.8 % 
transfers and 88.2 % recombinations). The ratio between the number of ruptures involved 
in a transfer and the total number of ruptures gives the recombination rate. 
Figs. IT5k and b give the percentage of the reactants and decay products for the O-H for- 
mations and ruptures including recombinations and Figs. Hot and d the related transfer 
reactions (i.e. without recombinations). For the transfers, due to the small numbers, the 
error bars are quite large. At both temperatures, the contributions of the Si-0 dangling 
bonds are higher for the transfer related ruptures than for the overall ruptures, i.e. transfers 
and recombinations. Furthermore we recognize from Figs. Hot and d that the Si-0 dan- 
gling bonds and the BOs are the main acceptors and rupture products, with almost equal 
probability. 

Concerning transfers, the Si-0 dangling bonds react with a probability of 40 % at 3500 K 
and 60 % at 3000 K into an SiOH group and the SiOH group (O*) decays with a probability 
of 40 % at 3500 K and 50 % at 3000 K into a Si-0 dangling bond. However the corresponding 
probability that includes recombinations is for both the formation and the decay around 30 
%. This result indicates that the Si-0 dangling bonds are more involved in the formation 
and rupture of the O-H bonds when the transfers are considered and the recombinations 
are taken out. Furthermore it also demonstrates the existence of large vibrations of the H 
atoms around the O atoms of the Si-0 dangling bonds and hence the existence of the "weak" 
hydrogen bonds described in Sec. IIII Al 

The BOs serve, with a probability of 50 % at 3500 K and of 40 % at 3000 K, as reactants 
for a hydrogen transfer leading to a new O-H group (60 % including recombinations at both 
temperatures). On the other side, the BOs serve as products of a hydrogen transfer with a 
probability of about 50 %, and therefore the participation of the BOs to the formation and 
rupture of the O-H bonds is slightly decreased when transfers are compared to the reactions 
including recombinations. 
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FIG. 16: Probabilities of O-H group formation and decay including recombinations, (a) and (b), 
and without recombinations, (c) and (d), associated to Si-0 dangling bonds, O* and BO. Black 
and white bars correspond to 3500 K and 3000 K, respectively. 



Finally, the O* contribute much less to the O-H group formation and decay than the other 
species, with a probability of about 10 % (12 % including recombinations) at 3500 K and 
% (8 % including recombinations) at 3000 K. 

From Fig. EE we deduce that hydrogen transfers involving Si-0 dangling bonds or BOs as 
reactants as well as decay products dominate with pa 90 %. Since these units are associated 
to O* or 03H as initial or final states, they correspond to the hydrogen diffusion process 
number 1. which means that among the processes 1. and 2., the former dominates with pa 
90 %. 

Let us now consider the third considered hydrogen diffusion mechanism. It is associated to 
the motion of a O-H unit as a whole (without any O-H rupture) which can be obtained via 
the formation and decay of an 03H cluster as sketched in Fig. El Note that this process 
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FIG. 17: Typical hydrogen diffusion reaction of process three involving a 03H. A O-H group 
is released from an Si-O-H (O*) to a 03H cluster decaying subsequently in another O*. An 
oversaturated and an unsaturated silicon atoms are subsequent products. 
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requires the presence of "weak" Si-0 bonds, the existence of which was inferred from Fig. 
El and leads to the formation of threefold and fivefold coordinated silicon atoms. As for 
the O-H bonds, we have analyzed the probabilities to find the different reactants and decay 
products for the 03H clusters formation and rupture. Along the equilibrated part of the 
trajectories, we have found 142 formations and 139 ruptures at 3500 K and 104 formations 
and 103 ruptures at 3000 K. In the case of 03H clusters, the possible reactants and decay 
products are the oxygen atoms involved in a Si-O-H unit (O*) or a bridging oxygen (BO). 
The latter case involves an O-H rupture event whereas the first case does not. This is the 
reason why a distinction between recombinations and transfers is not easily feasible for these 
species since the transfered object can either be a hydrogen atom (as in Tab. |Hj) or an entire 
O-H group (as in Fig. H7|) . We hence display the distribution for 03H formations and 
ruptures including recombinations in Fig. There are only very small differences between 
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FIG. 18: 03H group formation (left) and decay (right) products at 3500 K (black bars) and 3000 
K (white bars) including recombinations. 



formation and decay. At both temperatures, the 03H clusters form and decay from and into 
BO and O*, with almost the same probabilities. However at 3000 K, the probability from 
and into O* is slightly preferred with a value of 56 % for O* vs. 44 % for BO. This indicates 
that the hydrogen diffusion process involving the motion of a O-H group without any O-H 
bond rupture (process number 3.) occurs in the liquids as well. 

To conclude the section on the diffusion process, we analyzed the activation energy for a 
hydrogen transfer corresponding to the rupture of an O-H bond (without recombination). 
We can calculate the activation energy for a hyrogen transfer reaction assuming an Arrhenius 
behavior 

« (r) ^ex P {-^}, (6) 

where the reaction rate R(T) is given by a constant A and the activation energy Using 
the two numbers of O-H ruptures leading to transfers of 28 for 3500 K and 8 for 3000 K we ob- 
tain 280 kJ/mol (2.91 eV) for the activation energy (the different lengths of the equilibrated 
trajectories already taken into account). Unfortunately the constant A depends explicitly 
on the transfers per time and can hence not be determined due to the above mentioned 
problem of the thermostats. The activation energies for the generation and dissociation of 
various §i x O y }l z molecules in the gas phase have been determined by Zachariah et al |4j] 
with the GAUSSIAN ab initio package Zachariah et al found a value for the activation 
energy of the formation of OSiOH from Si02 and H of 285 kJ/mol which compares very well 
to the above given one of 280 kJ/mol for a hydrogen transfer. 
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IV. SUMMARY 



We have equilibrated a liquid of silica containing 3.84 wt. % H 2 at 3500 K and 3000 K 
using first-principles molecular-dynamics simulations. We observed that the water is mostly 
attached to the silica network in the form of Si-OH groups. Water molecules or free O-H 
groups occur only at the highest temperature but are not stable and decay rapidly. The 
SiC>4 tetrahedra are still the basic network forming units, as in pure silica. The short range 
correlation (i.e. the radial distribution functions) suggest that the structure of the matrix 
is as much changed by the addition of water than as by the addition of the same amount of 
sodium oxide to the liquid. In contrast to this, the way the modifier cation itself is attached 
to the matrix seems to be quite different. The sodium atoms tend to form bonds of ionic 
character with oxygen, whereas the hydrogen atoms are attached by covalent bonds. This 
difference in the bonding character of O-H and O-Na could be the reason for the slower 
diffusion of the hydrogen atoms in liquid silica compared to that of the sodium atoms. 
The neutron scattering structure factor shows a pronounced prepeak at a wave vector of 
about 1A _1 , i.e. we have evidence for the presence of long range correlations. The prepeak 
is even more intense than that for the sodium silicate at the same wave vector transfer. The 
origin of this feature in the hydrous liquid seems to be different from that of the sodium 
silicate, since the partial structure factor for H-H contributes only little. Furthermore, in 
the hydrous case, the silica matrix itself seems to be modified since the prepeak occurs as 
well in the partial structure factors for the matrix. It remains to be investigated whether the 
prepeak can be observed at ambient temperature. The data show a shift to higher g-vectors 
and a decrease of the intensity at the lowest temperature. 

The other experimentally accessible quantity, the Q n distribution, shows clearly that the 
O-H groups try to avoid each other. Silicon atoms with two and more Si-0 dangling bonds 
or O* (Q n sites with n < 2) are rarely found. Despite the strong covalent character of the 
O-H bond, the O-H units have a high radial mobility and oxygen triclusters have a high 
tendency to form. These oxygen triclusters violate the oxygen stoichiometry and the system 
compensates this violation by the formation of Si-0 dangling bonds. 

The dangling bonds and the triclusters constitute mainly the intermediate states for the 
diffusion mechanism of the hydrogen atoms in the melt. Indeed three possible diffusion 
mechanisms have been discussed, two of them involving the rupture of an O-H bond. Free 
water molecules, hydrogen molecules and free hydroxyl groups are not or only rarely found 
and cannot therefore be made responsible for the hydrogen diffusion in the liquid. 
Coming back to one of the motivating questions, the bubble formation, only little can be 
concluded from this study. On the one hand, it turned out that the O-H groups are very 
stable species and it requires a certain stoichiometry violation to weaken the O-H bonds 
(note that in O-H-0 transition units the stoichiometry of the H atom is also violated). On 
the other hand these bond weakening transition states occur in sufficiently high quantity 
to assure the hydrogen diffusion. Nevertheless, a clustering of hydrogen and the formation 
of stable water molecules do not take place, especially at the lower temperature. If a the 
formation of bubbles were realized, we would expect a higher signal in the H-H radial 
distribution function at typical H-H distances found in liquid water (2-3 A) and perhaps a 
signal in the Q n distribution for n < 2. Hence we do not see evidence for water clustering 
at the considered temperatures. We recognize that the influence of temperature, pressure 
and silicate composition on the bubble formation is not yet investigated and further work 
of the present kind will be mandatory. 
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